
%%%%this file uses permutations

global gammagrid elementsize lowerbound table m n agent  P beta S state agent Y scalefactor rho 
if beta<=0.93
elementsize=0.1;
lowerbound=0.1;
else
elementsize=0.01;
lowerbound=0.05;
end
[gammagrid,m]=grid3(elementsize,lowerbound);
agent=vss;

gammagrid(:,4)=1;
gammagrid(:,5)=gammagrid(:,2)./gammagrid(:,1);
gammagrid(:,6)=gammagrid(:,3)./gammagrid(:,1);
n=length(gammagrid);
table=[];
for r1=1:m
    index1=sum(m:-1:m-(r1-2));
    for r2=1:m-(r1-1)
        table(r1,r2)=index1+r2;
    end
end

d1=n;   %length of gammagrid
d2=8;  %states
d3=3;   %agent


WSB=ones(d1,d2,d3)*1000;
PHI=ones(d1,d2,d3)*1000;
C=ones(d1,d2,d3)*1000;

XXKEEP=ones(33,d2,d3)*1000;
GAMMAGRID=gammagrid;
PROB=ones(d2,d2)*1000;
STATE=ones(d2,d3)*1000;
TABLE=table;
COUNTGRID=ones(d1,d2);

income=incomevals;
S=gridmake(income,income,income);
SS=gridmake([1:length(income(:,1)) ]', [1:length(income(:,1)) ]', [1:length(income(:,1)) ]');

state=length(S);
Y=sum(S,2);				% possible aggregate income outcome

for k=1:state
    for j=1:state
        P(j,k)=Q1(SS(j,1),SS(k,1))*Q1(SS(j,2),SS(k,2))*Q1(SS(j,3),SS(k,3));
    end
end

% Compute set of autarky values
caut=zeros(state,agent);
for k=1:agent
    caut(:,k)=S(:,k);
end
waut=zeros(state,agent);
wautnew=zeros(state,agent),
dwaut=1;
crit=0.0000000001;
while dwaut>crit
    
    wautnew(:,1)=utility(rho,caut(:,1))+beta*P*waut(:,1);
    wautnew(:,2)=utility(rho,caut(:,2))+beta*P*waut(:,2);
    wautnew(:,3)=utility(rho,caut(:,3))+beta*P*waut(:,3);
    dwaut=max(max(abs(wautnew-waut)));
    waut=wautnew;
end


for j=1:state
    for g=1:agent
        if waut(j,g)>0
            waut(j,g)=waut(j,g)*(1-penalty);
        else
            waut(j,g)=waut(j,g)/(1-penalty);
        end
    end
end

% Solve problem without enforcement constraints at every grid point

w=zeros(n,state,agent);
gammar=zeros(n,state,agent);
c=zeros(n,state,agent);

% Note that phi contains the Lagrange multipliers on the enforcement
% constraint
phi=zeros(n,state,agent);

count=0;


for j=1:state
    for g=1:agent
        c(:,j,g)=Y(j)*gammagrid(:,agent+g)./(1+gammagrid(:,5)+gammagrid(:,6));
    end
end
cfirstbest=c;
gammar=zeros(n,state,agent);
for j=1:state
    for k=1:agent
        gammar(:,j,k)=gammagrid(:,agent+k);
        
    end
end
% Compute the value functions
w=zeros(n,state,agent);
wnew=zeros(n,state,agent);
dw=1;
while dw>crit
    for j=1:state
        wnew(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta*w(:,:,1)*P(j,:)';
        wnew(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta*w(:,:,2)*P(j,:)';
        wnew(:,j,3)=utility(rho,cfirstbest(:,j,3))+beta*w(:,:,3)*P(j,:)';

    end
            dw=max(max(max(abs(w-wnew))));
        w=wnew;
end
    wfirstbest=w;
    %%%
    ewfirstbest=zeros(n,state,agent);
    for j=1:state
        for g=1:agent
            ewfirstbest(:,j,g)=wfirstbest(:,:,g)*P(j,:)';
        end
    end
    ewsb=ewfirstbest;
    wsb=wfirstbest;
    
   
xkeep=[];
outputkeep=[];
xaux=[];
csol=[];
phisol=[];
wsbnewsol=[];
wsbnew=wsb;    
    gap=1;
    gapend=0.1;
    options=optimset('Display','iter');
    scalefactor=1;
    
    count=0;
    eps=0.0000001;
    t=0;
    gap=1;
    agentperm=perms(1:agent);
     gridperm=[]; iperm=[]; jperm=[]; csol=[]; wsbnewsol=[]; phisol=[]; 
    
    
     for t=1:25
        disp([gap,t])
        countgrid=ones(n,state);
        countgridod=ones(1,state);
        
        for j=1:state
            
            stateperm=perms(S(j,:));
            if countgridod(j)>0
            index=[1 j];
            clear x0
          
            %%for each state find the allocation if it binds for 1,2; 1,3; 2,3 so you
            %%don't need to repeat below:
            %%starting guesses
            rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<=waut(j,2));
            x0=[c(rr(end),j,1) c(rr(end),j,2)];
            x0=[x0 x0(2)/x0(1) (Y(j)-x0(2)-x0(1))/x0(1)];
            [xsol,output]=agent12solverho(x0,index,ewsb,waut);
            outputkeep12(j,1:3)=output(1:3);
            xkeep12(j,1:3)=[xsol(1:2) Y(j)-xsol(1)-xsol(2)];
            if output(6)<=0
                stop 
            end
            jjperm=find(S(:,1)==S(j,2) & S(:,2)==S(j,1) & S(:,3)==S(j,3));
            outputkeep12(jjperm,1)=outputkeep12(j,2);
            outputkeep12(jjperm,2)=outputkeep12(j,1);
            outputkeep12(jjperm,3)=outputkeep12(j,3);
            xkeep12(jjperm,1)=xkeep12(j,2);
            xkeep12(jjperm,2)=xkeep12(j,1);
            xkeep12(jjperm,3)=xkeep12(j,3);
            countgridod(jjperm)=0;
            end  
            jjjperm=find(S(:,1)==S(j,3) & S(:,2)==S(j,1) & S(:,3)==S(j,2));
            outputkeep23(jjjperm,1)=outputkeep12(j,3);
            outputkeep23(jjjperm,2)=outputkeep12(j,1);
            outputkeep23(jjjperm,3)=outputkeep12(j,2);

            xkeep23(jjjperm,1)=xkeep12(j,3);
            xkeep23(jjjperm,2)=xkeep12(j,1);
            xkeep23(jjjperm,3)=xkeep12(j,2);


            jjjperm=find(S(:,1)==S(j,1) & S(:,2)==S(j,3) & S(:,3)==S(j,2));
            outputkeep13(jjjperm,1)=outputkeep12(j,1);
            outputkeep13(jjjperm,2)=outputkeep12(j,3);
            outputkeep13(jjjperm,3)=outputkeep12(j,2);

            xkeep13(jjjperm,1)=xkeep12(j,1);
            xkeep13(jjjperm,2)=xkeep12(j,3);
            xkeep13(jjjperm,3)=xkeep12(j,2);

            
              
            
            
            
            
        end
        
        for j=1:state
        stateperm=perms(S(j,:));    
        for i=1:n
                
            index=[i,j];
            if countgrid(i,j)>0
            gridperm=perms(gammagrid(i,1:3));
            for k=1:length(gridperm)
            iperm(k)=find(gammagrid(:,1)==gridperm(k,1)& gammagrid(:,2)==gridperm(k,2)& gammagrid(:,3)==gridperm(k,3));
            jperm(k)=find(S(:,1)==stateperm(k,1) & S(:,2)==stateperm(k,2) & S(:,3)==stateperm(k,3));
            end
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%Constraint binding for  Person 1 only%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                if utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)<waut(j,1) ...
                        & utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ...
                        & utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)>=waut(j,3)
                    
                 
                    
                %%%%check whether it's possible that the individual deviation goes through    
                ccheck(2)=(Y(j)-xkeep12(j,1))/(1+gammagrid(i,3)/gammagrid(i,2));
                ccheck(3)=Y(j)-xkeep12(j,1)-ccheck(2);
                
                if ccheck(2)>=0.9*xkeep12(j,2) & ccheck(3)>=0.9*xkeep13(j,3)
                %%%auxiliary problem to get good starting guess
                clear x0
                xaux2 = (Y(j)-xkeep12(j,1))/(1+gammagrid(i,6)/gammagrid(i,5));
                xaux3 = gammagrid(i,6)/gammagrid(i,5)*xaux2;
                x0(1)=Y(j)-xaux2-xaux3;
                x0(2)=xaux2;
                x0=[x0 x0(2)/x0(1) (Y(j)-x0(2)-x0(1))/x0(1)];

                [xsol, output]=agent1solverho(x0,index,ewsb,waut);
                if output(6)<=0
                   error('did not converge') 
                end
                csol(i,j,1:2)=xsol(1:2);
                csol(i,j,3)=Y(j)-xsol(1)-xsol(2);
                gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                wsbnewsol(i,j,1:agent)=output(1:agent);
                phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
                phisol(i,j,2:agent)=0;
                %%%%%here check that it goes through ex-post
                if output(2)<waut(j,2) 
                
                csol(i,j,1)=xkeep12(j,1);
                csol(i,j,2)=xkeep12(j,2);
                phisol(i,j,3)=0;

                csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
                gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
                phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;
                wsbnewsol(i,j,:)=outputkeep12(j,1:agent);
                elseif output(3)<waut(j,3)
                csol(i,j,1:2)=xkeep13(j,1:2);
                csol(i,j,3)=Y(j)-xkeep13(j,1)-xkeep13(j,2);
                gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                wsbnewsol(i,j,1:agent)=outputkeep13(j,1:agent);
                phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
                phisol(i,j,3)=((gammar(i,j,3)/gammar(i,j,2))/(gammagrid(i,6)/gammagrid(i,5)))^rho-1;
                phisol(i,j,2)=0;    
                end
                %%%%here you check that it goes through ex-ante 
                elseif ccheck(2)<0.9*xkeep12(j,2) 
                
               
                csol(i,j,1)=xkeep12(j,1);
                csol(i,j,2)=xkeep12(j,2);
                phisol(i,j,3)=0;

                csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
                gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
                phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;
                wsbnewsol(i,j,:)=outputkeep12(j,1:agent);
                elseif ccheck(3)<0.9*xkeep13(j,3)
                
                csol(i,j,1:2)=xkeep13(j,1:2);
                csol(i,j,3)=Y(j)-xkeep13(j,1)-xkeep13(j,2);
                gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                wsbnewsol(i,j,1:agent)=outputkeep13(j,1:agent);
                phisol(i,j,1)=(gammagrid(i,5)/gammar(i,j,2))^rho-1;
                phisol(i,j,3)=((gammar(i,j,3)/gammar(i,j,2))/(gammagrid(i,6)/gammagrid(i,5)))^rho-1;
                phisol(i,j,2)=0;
               
                end
                    
                    %identical solutions
                    for k=1:length(gridperm)
                    for g=1:agent
                    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
                    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
                    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
                    end
                    countgrid(iperm(k),jperm(k))=0;
                    end
                    
                    
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Constraint binding for Person 1 and Person 2%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    elseif  utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)<waut(j,1) ...
                        & utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)<waut(j,2) ...
                        & utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)>=waut(j,3);
                    
                    
                    csol(i,j,1:2)=xkeep12(j,1:2);
                    csol(i,j,3)=Y(j)-csol(i,j,1)-csol(i,j,2);
                    phisol(i,j,3)=0;
                    
                    gammar(i,j,2)=csol(i,j,2)/csol(i,j,1);
                    gammar(i,j,3)=csol(i,j,3)/csol(i,j,1);
                    phisol(i,j,1)=(gammagrid(i,6)/gammar(i,j,3))^rho-1;
                    phisol(i,j,2)=((gammagrid(i,6)/gammagrid(i,5))/(gammar(i,j,3)/gammar(i,j,2)))^rho-1;
                    wsbnewsol(i,j,1:agent)=outputkeep12(j,1:agent);
                    
                       
                    %identical solutions
                    for k=1:length(gridperm)
                    for g=1:agent
                    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
                    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
                    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
                    end
                    countgrid(iperm(k),jperm(k))=0;
                    end
                    
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %Constraint binding for noone%%%%%%%
                    elseif utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)>=waut(j,1) ...
                     & utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ...
                     & utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3)>=waut(j,3);
                    csol(i,j,:)=cfirstbest(i,j,:);
                    gammar(i,j,:)=gammagrid(i,4:6);
                    phisol(i,j,1:agent)=0;
                    for g=1:agent
                        wsbnewsol(i,j,g)=utility(rho,cfirstbest(i,j,g))+beta*ewsb(i,j,g);
                    end
                 
                    %identical solutions
                    for k=1:length(gridperm)
                    for g=1:agent
                    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
                    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
                    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
                    end
                    countgrid(iperm(k),jperm(k))=0;
                    end
                end
                
                end %end of countgrid loop
            end % end of i loop
        end % end of j loop
        
        
        
        
        for g=1:agent
            gap(g)=scalefactor*norm((wsbnew(:,:,g)-wsb(:,:,g)));
        end
        for j=1:state
            for g=1:agent
                gammar(:,j,g)=c(:,j,g)./c(:,j,1);
            end
        end
        gap=sum(gap)
        [CC,II]=max(abs(wsbnew(:,:,1)-wsb(:,:,1)));
        wsb=wsbnew;
        
        for j=1:state
            for g=1:agent
                ewsb(:,j,g)=wsb(:,:,g)*P(j,:)';
            end
        end
        
    end % end of iteration loop
    
    
 
    
    
